General mixed linear models
Tracy Sánchez
Última actualización: 09 November 2022
Loading Data
X ID Alter Global Overview Cardinal AffiliativeQ_average item3
1 0 365 37 5.6 6.285714 2.5 5.2 2
2 1 3572 23 5.7 3.142857 1.0 4.0 2
3 2 3572 23 5.7 3.142857 1.0 4.0 2
4 3 3572 23 5.7 3.142857 1.0 4.0 2
5 4 3572 23 5.7 3.142857 1.0 4.0 2
6 5 3572 23 5.7 3.142857 1.0 4.0 2
ContextQ_average remember_action_likert remember_standing_likert
1 3.5 2.714286 1.107143
2 4.5 2.250000 1.714286
3 4.5 2.250000 1.714286
4 4.5 2.250000 1.714286
5 4.5 2.250000 1.714286
6 4.5 2.250000 1.714286
attractive_action_likert attractive_standing_likert realistic_action_likert
1 2.214286 1.178571 3.392857
2 2.607143 2.607143 2.964286
3 2.607143 2.607143 2.964286
4 2.607143 2.607143 2.964286
5 2.607143 2.607143 2.964286
6 2.607143 2.607143 2.964286
realistic_standing_likert AbsolutError SignedAngle..180 IQR RT
1 1.964286 NA NA NA NA
2 2.678571 34.308001 -34.308001 52.35318 16.368489
3 2.678571 83.069424 83.069424 52.35318 7.974188
4 2.678571 2.210407 2.210407 52.35318 8.858059
5 2.678571 30.482306 -30.482306 52.35318 6.429150
6 2.678571 5.407388 5.407388 52.35318 4.768045
DistanceToParticipant PointingTaskStartingLocations TrialNumber StartPointID
1 NA NA NA NA
2 263.1861 18 19 2
3 297.9136 7 274 23
4 173.0462 22 282 24
5 228.5974 12 293 25
6 273.0915 9 330 28
ID_for_StartingPosition avatarID AvatarPresenceCategory Context
1 NA NA
2 7 5 Omitted True
3 10 5 Present True
4 6 5 Present True
5 5 5 Present True
6 6 5 Omitted True
meaningfulBuilding ImageName
1
2 Meaningful 05_CmANo
3 Meaningful 05_CmA
4 Meaningful 05_CmA
5 Meaningful 05_CmA
6 Meaningful 05_CmANoRecoding for effects and formating factors
HumanA$ContextEffectf <-dplyr::recode(HumanA$Context,
False = -0.5, True= 0.5,
.default = NaN)
HumanA$AgentPresence <-dplyr::recode(HumanA$AvatarPresenceCategory,
Omitted = -0.5, Present= 0.5,
.default = NaN)
HumanA$ContextEffectf <-factor(HumanA$ContextEffectf,levels= c(-0.5, 0.5),
labels=c('Residential', 'Public'))
HumanA$AgentPresencef <-factor(HumanA$AgentPresence,
levels= c(-0.5, 0.5),
labels=c('Omitted', 'Displayed')) Data descriptives
MainVariables <- subset(HumanA, select = c(AbsolutError, RT))
summary(MainVariables) AbsolutError RT
Min. : 0.00986 Min. : 1.061
1st Qu.: 12.43032 1st Qu.: 3.845
Median : 33.48641 Median : 6.438
Mean : 48.17395 Mean : 8.080
3rd Qu.: 71.75575 3rd Qu.:10.717
Max. :179.98339 Max. :29.366
NA's :4 NA's :4 summary(HumanA$AgentPresence) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-0.500000 -0.500000 -0.500000 -0.000426 0.500000 0.500000 4 df = HumanA[complete.cases(HumanA),]
ggplot(df, aes(x=ContextEffectf, y=AbsolutError, fill=AgentPresencef)) +
geom_boxplot(notch=TRUE,
notchwidth = 0.8,
outlier.colour="red",
outlier.fill="red",
outlier.size=0.5)
ggplot(df, aes(x=ContextEffectf, y=RT, fill=AgentPresencef)) +
geom_boxplot(notch=TRUE,
notchwidth = 0.8,
outlier.colour="red",
outlier.fill="red",
outlier.size=0.5)
library(dplyr)
TwoFactorTable <- HumanA %>%
group_by(ContextEffectf, AgentPresencef)%>%
summarise(AccuracyMean = mean(AbsolutError, na.rm = TRUE),
n=n(),
AccuracyStandardDev = sd(AbsolutError, na.rm = TRUE),
RTMean = mean(RT, na.rm = TRUE),
RTStandardDev = sd(RT, na.rm = TRUE))`summarise()` has grouped output by 'ContextEffectf'. You can override using
the `.groups` argument.library(tidyr)
Attaching package: 'tidyr'
The following objects are masked from 'package:Matrix':
expand, pack, unpack
The following object is masked from 'package:dlookr':
extractTwoFactorTableUnite <- TwoFactorTable %>%
unite("TwoFactor", ContextEffectf:AgentPresencef, sep= " ", remove = F)
TwoFactorTableUnite <- TwoFactorTableUnite %>%
mutate( AccuracyStandardError=AccuracyStandardDev/sqrt(n)) %>%
mutate( AccuracyStandardIC=AccuracyStandardDev * qt((1-0.05)/2 + .5, n-1)) %>%
mutate( RTStandardError=RTStandardDev/sqrt(n)) Warning: Ignoring unknown aesthetics: linetype
Warning: Ignoring unknown aesthetics: linetype
Checking for the distribution of Absolut error
df = HumanA[complete.cases(HumanA),]
df$AbsolutErrorR <- round(df$AbsolutError, digits = 3)
qqp(df$AbsolutErrorR, "norm")
[1] 2613 968qqp(df$AbsolutErrorR, "lnorm")
[1] 2613 968Linear mixed models
Assessing the need for a multilevel model
interceptOnly <-gls(AbsolutError ~ 1, data = df,
method = "ML")
IDrandomInterceptOnly <-lme(AbsolutError ~ 1, data = df,
random =~1|ID,
method = "ML")
StartlocationsrandomIntercept <-update(IDrandomInterceptOnly, .~.,
random=~1|ID/PointingTaskStartingLocations,
method= "ML")Including Id and starting position as random effects significantly improves the fit of the model
Absolut Error Models
I am adding one main factor at a time
MeaningfulContext <-update(StartlocationsrandomIntercept, .~. + ContextEffectf)
summary(MeaningfulContext)Linear mixed-effects model fit by maximum likelihood
Data: df
AIC BIC logLik
90086.83 90122.2 -45038.42
Random effects:
Formula: ~1 | ID
(Intercept)
StdDev: 13.63238
Formula: ~1 | PointingTaskStartingLocations %in% ID
(Intercept) Residual
StdDev: 13.67445 40.87166
Fixed effects: AbsolutError ~ ContextEffectf
Value Std.Error DF t-value p-value
(Intercept) 50.85647 2.7933018 7984 18.206580 0.0000
ContextEffectfPublic -2.51541 0.9064463 7984 -2.775021 0.0055
Correlation:
(Intr)
ContextEffectfPublic -0.162
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.1474521 -0.6819886 -0.2313140 0.5171520 3.5682778
Number of Observations: 8713
Number of Groups:
ID PointingTaskStartingLocations %in% ID
26 728 Anova(MeaningfulContext)Analysis of Deviance Table (Type II tests)
Response: AbsolutError
Chisq Df Pr(>Chisq)
ContextEffectf 7.7025 1 0.005514 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Presence <-update(MeaningfulContext, .~. + AgentPresencef)
summary(Presence)Linear mixed-effects model fit by maximum likelihood
Data: df
AIC BIC logLik
90087.93 90130.37 -45037.97
Random effects:
Formula: ~1 | ID
(Intercept)
StdDev: 13.63279
Formula: ~1 | PointingTaskStartingLocations %in% ID
(Intercept) Residual
StdDev: 13.66899 40.87023
Fixed effects: AbsolutError ~ ContextEffectf + AgentPresencef
Value Std.Error DF t-value p-value
(Intercept) 50.43022 2.8294125 7983 17.823567 0.0000
ContextEffectfPublic -2.51960 0.9064670 7983 -2.779579 0.0055
AgentPresencefDisplayed 0.85750 0.9039836 7983 0.948584 0.3429
Correlation:
(Intr) CntxEP
ContextEffectfPublic -0.159
AgentPresencefDisplayed -0.159 -0.005
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.1528440 -0.6787019 -0.2308530 0.5148011 3.5789589
Number of Observations: 8713
Number of Groups:
ID PointingTaskStartingLocations %in% ID
26 728 Anova(Presence)Analysis of Deviance Table (Type II tests)
Response: AbsolutError
Chisq Df Pr(>Chisq)
ContextEffectf 7.7287 1 0.005435 **
AgentPresencef 0.9001 1 0.342749
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TwofactorInteraction <-update(Presence, .~. + ContextEffectf*AgentPresencef)
summary(TwofactorInteraction)Linear mixed-effects model fit by maximum likelihood
Data: df
AIC BIC logLik
90089.62 90139.13 -45037.81
Random effects:
Formula: ~1 | ID
(Intercept)
StdDev: 13.62648
Formula: ~1 | PointingTaskStartingLocations %in% ID
(Intercept) Residual
StdDev: 13.67256 40.86898
Fixed effects: AbsolutError ~ ContextEffectf + AgentPresencef + ContextEffectf:AgentPresencef
Value Std.Error DF t-value
(Intercept) 50.18512 2.861825 7982 17.536055
ContextEffectfPublic -2.01680 1.273328 7982 -1.583881
AgentPresencefDisplayed 1.35775 1.268506 7982 1.070352
ContextEffectfPublic:AgentPresencefDisplayed -1.01262 1.801138 7982 -0.562209
p-value
(Intercept) 0.0000
ContextEffectfPublic 0.1133
AgentPresencefDisplayed 0.2845
ContextEffectfPublic:AgentPresencefDisplayed 0.5740
Correlation:
(Intr) CntxEP AgntPD
ContextEffectfPublic -0.219
AgentPresencefDisplayed -0.219 0.490
ContextEffectfPublic:AgentPresencefDisplayed 0.152 -0.702 -0.702
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.1599911 -0.6810745 -0.2312981 0.5155331 3.5727542
Number of Observations: 8713
Number of Groups:
ID PointingTaskStartingLocations %in% ID
26 728 Anova(TwofactorInteraction)Analysis of Deviance Table (Type II tests)
Response: AbsolutError
Chisq Df Pr(>Chisq)
ContextEffectf 7.7287 1 0.005435 **
AgentPresencef 0.9001 1 0.342765
ContextEffectf:AgentPresencef 0.3162 1 0.573886
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pairwise<- emmeans(TwofactorInteraction, pairwise ~ ContextEffectf*AgentPresencef)
Pairwise$emmeans
ContextEffectf AgentPresencef emmean SE df lower.CL upper.CL
Residential Omitted 50.2 2.86 25 44.3 56.1
Public Omitted 48.2 2.87 25 42.3 54.1
Residential Displayed 51.5 2.87 25 45.6 57.4
Public Displayed 48.5 2.86 25 42.6 54.4
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Residential Omitted - Public Omitted 2.017 1.27 7982 1.584 0.3879
Residential Omitted - Residential Displayed -1.358 1.27 7982 -1.070 0.7077
Residential Omitted - Public Displayed 1.672 1.28 7982 1.309 0.5573
Public Omitted - Residential Displayed -3.375 1.28 7982 -2.630 0.0426
Public Omitted - Public Displayed -0.345 1.28 7982 -0.269 0.9932
Residential Displayed - Public Displayed 3.029 1.28 7982 2.363 0.0845
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 4 estimates plot(Pairwise[[2]], CIs = TRUE)
library(multcomp);library(multcompView)Loading required package: mvtnormLoading required package: survivalLoading required package: TH.data
Attaching package: 'TH.data'The following object is masked from 'package:MASS':
geyserCLD <- cld(Pairwise,
alpha=0.05,
Letters=letters,
adjust="sidak")I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
See '? emm_list' for more informationggplot(CLD,
aes(x = ContextEffectf,
y = emmean,
group = AgentPresencef,
colours = .group)) +
geom_point(aes(shape=AgentPresencef, linetype =AgentPresencef), position=position_dodge(0.3)) +
geom_errorbar(aes(linetype=AgentPresencef,
ymin = lower.CL,
ymax = upper.CL),
position=position_dodge(0.3),
width = 0.2,
size = 0.7) +
theme_bw() +
theme(axis.title = element_text(face = "bold"),
axis.text = element_text(face = "bold"),
plot.caption = element_text(hjust = 0)) +
ylab("Estimated marginal mean\ Absolute angular error") +
xlab("Location location") +
ggtitle ("Marginal Means",
subtitle = "location * Presence") +
labs(caption = paste0(
"Boxes indicate the EM mean. \n",
"Error bars indicate the 95% ",
"confidence interval of the EM mean. \n"),
hjust=0.5) Warning: Ignoring unknown aesthetics: linetype
anova(interceptOnly, IDrandomInterceptOnly, StartlocationsrandomIntercept,
MeaningfulContext, Presence, TwofactorInteraction ) Model df AIC BIC logLik Test
interceptOnly 1 2 91151.27 91165.42 -45573.64
IDrandomInterceptOnly 2 3 90386.56 90407.77 -45190.28 1 vs 2
StartlocationsrandomIntercept 3 4 90092.53 90120.82 -45042.27 2 vs 3
MeaningfulContext 4 5 90086.83 90122.20 -45038.42 3 vs 4
Presence 5 6 90087.93 90130.37 -45037.97 4 vs 5
TwofactorInteraction 6 7 90089.62 90139.13 -45037.81 5 vs 6
L.Ratio p-value
interceptOnly
IDrandomInterceptOnly 766.7188 <.0001
StartlocationsrandomIntercept 296.0217 <.0001
MeaningfulContext 7.6986 0.0055
Presence 0.9000 0.3428
TwofactorInteraction 0.3162 0.5739plot(TwofactorInteraction, which = 1)
plot(MeaningfulContext, which = 1)
GHQ <- glmer(AbsolutError ~ ContextEffectf*AgentPresencef + (1|ID), data = HumanA,family=gaussian(link = "log"), nAGQ = 25) Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?summary(GHQ)Generalized linear mixed model fit by maximum likelihood (Adaptive
Gauss-Hermite Quadrature, nAGQ = 25) [glmerMod]
Family: gaussian ( log )
Formula: AbsolutError ~ ContextEffectf * AgentPresencef + (1 | ID)
Data: HumanA
AIC BIC logLik deviance df.resid
16766820 16766863 -8383404 16766808 9378
Scaled residuals:
Min 1Q Median 3Q Max
-1.6944 -0.7128 -0.2587 0.4884 3.7236
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 238.1 15.43
Residual 1786.7 42.27
Number of obs: 9384, groups: ID, 28
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.8286440 0.0690074 55.482
ContextEffectfPublic -0.0489551 0.0005874 -83.348
AgentPresencefDisplayed 0.0209055 0.0005661 36.929
ContextEffectfPublic:AgentPresencefDisplayed 0.0050108 0.0008215 6.099
Pr(>|z|)
(Intercept) < 2e-16 ***
ContextEffectfPublic < 2e-16 ***
AgentPresencefDisplayed < 2e-16 ***
ContextEffectfPublic:AgentPresencefDisplayed 1.06e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) CntxEP AgntPD
CntxtEffctP -0.004
AgntPrsncfD -0.004 0.490
CntxtEP:APD 0.003 -0.715 -0.690
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?Anova(GHQ)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: AbsolutError
Chisq Df Pr(>Chisq)
ContextEffectf 12764.563 1 < 2.2e-16 ***
AgentPresencef 3227.441 1 < 2.2e-16 ***
ContextEffectf:AgentPresencef 37.203 1 1.064e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(GHQ)Analysis of Variance Table
npar Sum Sq Mean Sq F value
ContextEffectf 1 12731.2 12731.2 7.1255
AgentPresencef 1 3234.7 3234.7 1.8104
ContextEffectf:AgentPresencef 1 37.2 37.2 0.0208emmeans(GHQ, pairwise ~ ContextEffectf:AgentPresencef, type = "response")$emmeans
ContextEffectf AgentPresencef response SE df asymp.LCL asymp.UCL
Residential Omitted 46.0 3.17 Inf 40.2 52.7
Public Omitted 43.8 3.02 Inf 38.3 50.1
Residential Displayed 47.0 3.24 Inf 41.0 53.8
Public Displayed 45.0 3.10 Inf 39.3 51.5
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null z.ratio
Residential Omitted / Public Omitted 1.0502 0.0006168 Inf 1 83.348
Residential Omitted / Residential Displayed 0.9793 0.0005544 Inf 1 -36.929
Residential Omitted / Public Displayed 1.0233 0.0005919 Inf 1 39.831
Public Omitted / Residential Displayed 0.9325 0.0005435 Inf 1 -119.859
Public Omitted / Public Displayed 0.9744 0.0005796 Inf 1 -43.567
Residential Displayed / Public Displayed 1.0449 0.0006001 Inf 1 76.513
p.value
<.0001
<.0001
<.0001
<.0001
<.0001
<.0001
P value adjustment: tukey method for comparing a family of 4 estimates
Tests are performed on the log scale plot(fitted(GHQ), residuals(GHQ), xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, lty = 2)
lines(smooth.spline(fitted(GHQ), residuals(GHQ)))
interceptOnly <-gls(log(AbsolutError) ~ 1, data = df,
method = "ML")
IDrandomInterceptOnly <-lme(log(AbsolutError) ~ 1, data = df,
random =~1|ID,
method = "ML")
StartlocationsrandomIntercept <-update(IDrandomInterceptOnly, .~.,
random=~1|ID/PointingTaskStartingLocations,
method= "ML")
MeaningfulContext <-update(StartlocationsrandomIntercept, .~. + ContextEffectf)
Presence <-update(MeaningfulContext, .~. + AgentPresencef)
TwofactorInteraction <-update(Presence, .~. + ContextEffectf*AgentPresencef)
summary(TwofactorInteraction)Linear mixed-effects model fit by maximum likelihood
Data: df
AIC BIC logLik
29477.89 29527.39 -14731.94
Random effects:
Formula: ~1 | ID
(Intercept)
StdDev: 0.4164386
Formula: ~1 | PointingTaskStartingLocations %in% ID
(Intercept) Residual
StdDev: 0.3017131 1.27877
Fixed effects: log(AbsolutError) ~ ContextEffectf + AgentPresencef + ContextEffectf:AgentPresencef
Value Std.Error DF t-value
(Intercept) 3.312583 0.08694747 7982 38.09867
ContextEffectfPublic -0.087901 0.03951064 7982 -2.22473
AgentPresencefDisplayed 0.064008 0.03939876 7982 1.62461
ContextEffectfPublic:AgentPresencefDisplayed -0.045421 0.05590503 7982 -0.81247
p-value
(Intercept) 0.0000
ContextEffectfPublic 0.0261
AgentPresencefDisplayed 0.1043
ContextEffectfPublic:AgentPresencefDisplayed 0.4165
Correlation:
(Intr) CntxEP AgntPD
ContextEffectfPublic -0.224
AgentPresencefDisplayed -0.224 0.491
ContextEffectfPublic:AgentPresencefDisplayed 0.156 -0.704 -0.703
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-5.8952345 -0.4486963 0.1935545 0.6817024 2.2375050
Number of Observations: 8713
Number of Groups:
ID PointingTaskStartingLocations %in% ID
26 728 Anova(TwofactorInteraction)Analysis of Deviance Table (Type II tests)
Response: log(AbsolutError)
Chisq Df Pr(>Chisq)
ContextEffectf 15.5081 1 8.215e-05 ***
AgentPresencef 2.1956 1 0.1384
ContextEffectf:AgentPresencef 0.6604 1 0.4164
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pairwise<- emmeans(TwofactorInteraction, pairwise ~ ContextEffectf*AgentPresencef, type='response')
Pairwise $emmeans
ContextEffectf AgentPresencef response SE df lower.CL upper.CL
Residential Omitted 27.5 2.39 25 23.0 32.8
Public Omitted 25.1 2.19 25 21.0 30.1
Residential Displayed 29.3 2.55 25 24.5 35.0
Public Displayed 25.6 2.23 25 21.4 30.6
Degrees-of-freedom method: containment
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio
Residential Omitted / Public Omitted 1.092 0.0431 7982 1 2.225
Residential Omitted / Residential Displayed 0.938 0.0370 7982 1 -1.625
Residential Omitted / Public Displayed 1.072 0.0423 7982 1 1.754
Public Omitted / Residential Displayed 0.859 0.0342 7982 1 -3.817
Public Omitted / Public Displayed 0.982 0.0390 7982 1 -0.468
Residential Displayed / Public Displayed 1.143 0.0454 7982 1 3.357
p.value
0.1166
0.3647
0.2957
0.0008
0.9661
0.0044
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 4 estimates
Tests are performed on the log scale plot(Pairwise[[2]])
anova(interceptOnly, IDrandomInterceptOnly, StartlocationsrandomIntercept,
MeaningfulContext, Presence, TwofactorInteraction ) Model df AIC BIC logLik Test
interceptOnly 1 2 30340.42 30354.56 -15168.21
IDrandomInterceptOnly 2 3 29586.02 29607.24 -14790.01 1 vs 2
StartlocationsrandomIntercept 3 4 29490.14 29518.44 -14741.07 2 vs 3
MeaningfulContext 4 5 29476.74 29512.10 -14733.37 3 vs 4
Presence 5 6 29476.55 29518.98 -14732.27 4 vs 5
TwofactorInteraction 6 7 29477.89 29527.39 -14731.94 5 vs 6
L.Ratio p-value
interceptOnly
IDrandomInterceptOnly 756.3928 <.0001
StartlocationsrandomIntercept 97.8778 <.0001
MeaningfulContext 15.4029 0.0001
Presence 2.1953 0.1384
TwofactorInteraction 0.6602 0.4165plot(TwofactorInteraction, which = 1)
model <- glmmPQL(RT ~ ContextEffectf*AgentPresencef, ~1|ID/PointingTaskStartingLocations, family = gaussian(link = "log"), data = HumanA, verbose = FALSE)
summary(model)Linear mixed-effects model fit by maximum likelihood
Data: HumanA
AIC BIC logLik
NA NA NA
Random effects:
Formula: ~1 | ID
(Intercept)
StdDev: 0.2231925
Formula: ~1 | PointingTaskStartingLocations %in% ID
(Intercept) Residual
StdDev: 0.2224999 4.923456
Variance function:
Structure: fixed weights
Formula: ~invwt
Fixed effects: RT ~ ContextEffectf * AgentPresencef
Value Std.Error DF
(Intercept) 2.0659373 0.04469080 8597
ContextEffectfPublic 0.0235019 0.01708111 8597
AgentPresencefDisplayed -0.0110505 0.01718929 8597
ContextEffectfPublic:AgentPresencefDisplayed -0.0217719 0.02441735 8597
t-value p-value
(Intercept) 46.22735 0.0000
ContextEffectfPublic 1.37590 0.1689
AgentPresencefDisplayed -0.64287 0.5203
ContextEffectfPublic:AgentPresencefDisplayed -0.89166 0.3726
Correlation:
(Intr) CntxEP AgntPD
ContextEffectfPublic -0.190
AgentPresencefDisplayed -0.188 0.489
ContextEffectfPublic:AgentPresencefDisplayed 0.131 -0.698 -0.700
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.2369429 -0.6791802 -0.3101101 0.3731269 4.3832247
Number of Observations: 9384
Number of Groups:
ID PointingTaskStartingLocations %in% ID
28 784 Response Time models
Checking for the distribution of RT
df = HumanA[complete.cases(HumanA),]
df$RTr <- round(df$RT, digits = 3)
qqp(df$RT, "norm")
[1] 8633 8491qqp(df$RT, "lnorm")
[1] 8633 8491interceptOnlyt <-gls(log(RT) ~ 1, data = df,
method = "ML")
IDrandomInterceptOnlyt <-lme(log(RT) ~ 1, data = df,
random =~1|ID,
method = "ML")
StartlocationsrandomInterceptt <-lme(log(RT) ~ 1, data = df,
random=~1|ID|PointingTaskStartingLocations,
method= "ML")
MeaningfulContext <-update(StartlocationsrandomInterceptt, .~. + ContextEffectf)
Presence <-update(MeaningfulContext, .~. + AgentPresencef)
TwofactorInteraction <-update(Presence, .~. + ContextEffectf*AgentPresencef)
summary(TwofactorInteraction)Linear mixed-effects model fit by maximum likelihood
Data: df
AIC BIC logLik
18247.11 18303.69 -9115.555
Random effects:
Formula: ~1 | ID | PointingTaskStartingLocations
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.08246629 (Intr)
1 | IDTRUE 0.08246629 -0.587
Residual 0.68713429
Fixed effects: log(RT) ~ ContextEffectf + AgentPresencef + ContextEffectf:AgentPresencef
Value Std.Error DF
(Intercept) 1.8453072 0.02037131 8682
ContextEffectfPublic 0.0364076 0.02086710 8682
AgentPresencefDisplayed -0.0228840 0.02082944 8682
ContextEffectfPublic:AgentPresencefDisplayed -0.0168581 0.02949335 8682
t-value p-value
(Intercept) 90.58364 0.0000
ContextEffectfPublic 1.74474 0.0811
AgentPresencefDisplayed -1.09864 0.2720
ContextEffectfPublic:AgentPresencefDisplayed -0.57159 0.5676
Correlation:
(Intr) CntxEP AgntPD
ContextEffectfPublic -0.504
AgentPresencefDisplayed -0.504 0.493
ContextEffectfPublic:AgentPresencefDisplayed 0.356 -0.706 -0.706
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.689589261 -0.748955546 -0.002091011 0.739456135 2.353193766
Number of Observations: 8713
Number of Groups: 28 Anova(TwofactorInteraction)Analysis of Deviance Table (Type II tests)
Response: log(RT)
Chisq Df Pr(>Chisq)
ContextEffectf 3.5895 1 0.05815 .
AgentPresencef 4.5065 1 0.03377 *
ContextEffectf:AgentPresencef 0.3269 1 0.56751
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pairwiset<- emmeans(TwofactorInteraction, pairwise ~ ContextEffectf*AgentPresencef, type='response')
Pairwiset $emmeans
ContextEffectf AgentPresencef response SE df lower.CL upper.CL
Residential Omitted 6.33 0.129 27 6.07 6.60
Public Omitted 6.56 0.135 27 6.29 6.85
Residential Displayed 6.19 0.127 27 5.93 6.45
Public Displayed 6.31 0.129 27 6.05 6.58
Degrees-of-freedom method: containment
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio
Residential Omitted / Public Omitted 0.964 0.0201 8682 1 -1.745
Residential Omitted / Residential Displayed 1.023 0.0213 8682 1 1.099
Residential Omitted / Public Displayed 1.003 0.0208 8682 1 0.161
Public Omitted / Residential Displayed 1.061 0.0223 8682 1 2.823
Public Omitted / Public Displayed 1.041 0.0217 8682 1 1.904
Residential Displayed / Public Displayed 0.981 0.0205 8682 1 -0.936
p.value
0.3006
0.6904
0.9985
0.0246
0.2265
0.7854
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 4 estimates
Tests are performed on the log scale ref_grid(TwofactorInteraction)'emmGrid' object with variables:
ContextEffectf = Residential, Public
AgentPresencef = Omitted, Displayed
Transformation: "log" plot(Pairwiset[[2]])
anova(interceptOnlyt, IDrandomInterceptOnlyt, StartlocationsrandomInterceptt,
MeaningfulContext, Presence, TwofactorInteraction ) Model df AIC BIC logLik Test
interceptOnlyt 1 2 18303.41 18317.55 -9149.703
IDrandomInterceptOnlyt 2 3 16981.94 17003.16 -8487.971 1 vs 2
StartlocationsrandomInterceptt 3 5 18249.43 18284.79 -9119.716 2 vs 3
MeaningfulContext 4 6 18247.94 18290.38 -9117.971 3 vs 4
Presence 5 7 18245.44 18294.94 -9115.719 4 vs 5
TwofactorInteraction 6 8 18247.11 18303.69 -9115.555 5 vs 6
L.Ratio p-value
interceptOnlyt
IDrandomInterceptOnlyt 1323.4635 <.0001
StartlocationsrandomInterceptt 1263.4883 <.0001
MeaningfulContext 3.4887 0.0618
Presence 4.5052 0.0338
TwofactorInteraction 0.3268 0.5676plot(TwofactorInteraction, which = 1)
CLD <- cld(Pairwiset,
alpha=0.05,
Letters=letters,
adjust="sidak")I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
See '? emm_list' for more informationggplot(CLD,
aes(x = ContextEffectf,
y = response,
group = AgentPresencef,
colours = .group)) +
geom_point(aes(shape=AgentPresencef, color=AgentPresencef), position=position_dodge(0.3)) +
geom_errorbar(aes(color=AgentPresencef,
ymin = lower.CL,
ymax = upper.CL),
position=position_dodge(0.3),
width = 0.2,
size = 0.7) +
theme_bw() +
theme(axis.title = element_text(face = "bold"),
axis.text = element_text(face = "bold"),
plot.caption = element_text(hjust = 0)) +
ylab("Estimated marginal mean\ Response Time in Seconds") +
xlab("Location location") +
ggtitle ("Marginal Means",
subtitle = "location * Presence") +
labs(caption = paste0(
"Boxes indicate the EM mean. \n",
"Error bars indicate the 95% ",
"confidence interval of the EM mean. \n"),
hjust=0.5) 
GHQ <- glmer(RT ~ ContextEffectf*AgentPresencef + (1|ID), data = HumanA,
family=gaussian(link = "log"), nAGQ = 25)
summary(GHQ)Generalized linear mixed model fit by maximum likelihood (Adaptive
Gauss-Hermite Quadrature, nAGQ = 25) [glmerMod]
Family: gaussian ( log )
Formula: RT ~ ContextEffectf * AgentPresencef + (1 | ID)
Data: HumanA
AIC BIC logLik deviance df.resid
262057.8 262100.7 -131022.9 262045.8 9378
Scaled residuals:
Min 1Q Median 3Q Max
-1.9259 -0.6951 -0.2851 0.4426 4.0698
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 1.581 1.257
Residual 27.904 5.282
Number of obs: 9384, groups: ID, 28
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.0612428 0.0450587 45.746
ContextEffectfPublic 0.0182153 0.0034898 5.220
AgentPresencefDisplayed -0.0009493 0.0035211 -0.270
ContextEffectfPublic:AgentPresencefDisplayed -0.0337878 0.0049858 -6.777
Pr(>|z|)
(Intercept) < 2e-16 ***
ContextEffectfPublic 1.79e-07 ***
AgentPresencefDisplayed 0.787
ContextEffectfPublic:AgentPresencefDisplayed 1.23e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) CntxEP AgntPD
CntxtEffctP -0.038
AgntPrsncfD -0.038 0.491
CntxtEP:APD 0.027 -0.701 -0.707Anova(GHQ)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: RT
Chisq Df Pr(>Chisq)
ContextEffectf 0.4343 1 0.5099
AgentPresencef 51.3019 1 7.920e-13 ***
ContextEffectf:AgentPresencef 45.9253 1 1.228e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(GHQ)Analysis of Variance Table
npar Sum Sq Mean Sq F value
ContextEffectf 1 0.359 0.359 0.0129
AgentPresencef 1 51.466 51.466 1.8444
ContextEffectf:AgentPresencef 1 46.000 46.000 1.6485emmeans(GHQ, pairwise ~ ContextEffectf:AgentPresencef )$emmeans
ContextEffectf AgentPresencef emmean SE df asymp.LCL asymp.UCL
Residential Omitted 2.06 0.0451 Inf 1.97 2.15
Public Omitted 2.08 0.0451 Inf 1.99 2.17
Residential Displayed 2.06 0.0451 Inf 1.97 2.15
Public Displayed 2.04 0.0451 Inf 1.96 2.13
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio
Residential Omitted - Public Omitted -0.018215 0.00349 Inf -5.220
Residential Omitted - Residential Displayed 0.000949 0.00352 Inf 0.270
Residential Omitted - Public Displayed 0.016522 0.00350 Inf 4.715
Public Omitted - Residential Displayed 0.019165 0.00354 Inf 5.419
Public Omitted - Public Displayed 0.034737 0.00352 Inf 9.856
Residential Displayed - Public Displayed 0.015573 0.00356 Inf 4.379
p.value
<.0001
0.9932
<.0001
<.0001
<.0001
0.0001
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 4 estimates HumanA$Agent_Category <- with(HumanA, ave(seq_along(ID), ID, FUN = function(x) sample(c(rep('Action', ceiling(length(x)*0.6)), rep('Standing', length(x) - ceiling(length(x)*0.6))))))library(dplyr)
TwoFactorTable <- HumanA %>%
group_by(ContextEffectf, AgentPresencef, Agent_Category)%>%
summarise(AccuracyMean = mean(AbsolutError, na.rm = TRUE),
n=n(),
AccuracyStandardDev = sd(AbsolutError, na.rm = TRUE),
RTMean = mean(RT, na.rm = TRUE),
RTStandardDev = sd(RT, na.rm = TRUE))`summarise()` has grouped output by 'ContextEffectf', 'AgentPresencef'. You can
override using the `.groups` argument.library(tidyr)
TwoFactorTableUnite <- TwoFactorTable %>%
unite("TwoFactor", ContextEffectf:Agent_Category, sep= " ", remove = F)
TwoFactorTableUnite <- TwoFactorTableUnite %>%
mutate( AccuracyStandardError=AccuracyStandardDev/sqrt(n)) %>%
mutate( AccuracyStandardIC=AccuracyStandardDev * qt((1-0.05)/2 + .5, n-1)) %>%
mutate( RTStandardError=RTStandardDev/sqrt(n)) ggplot(data=subset(TwoFactorTableUnite, !is.na(AgentPresencef)),
aes(x = ContextEffectf,
y = AccuracyMean,
group = Agent_Category)) +
geom_point(aes(shape=Agent_Category, linetype =Agent_Category), position=position_dodge(0.3)) +
geom_errorbar(aes(linetype=Agent_Category,
ymin= AccuracyMean-AccuracyStandardError,
ymax=AccuracyMean+AccuracyStandardError),
position=position_dodge(0.3),
width = 0.2,
size = 0.7) +
facet_wrap(~ AgentPresencef) +
theme_bw() +
theme(axis.title = element_text(face = "bold"),
axis.text = element_text(face = "bold"),
plot.caption = element_text(hjust = 0)) +
ylab("Absolute angular error in degrees") +
xlab("Location location") +
ggtitle ("Accuracy performance",
subtitle = "The effect of location and presence") +
labs(caption = paste0(
"Error bars indicate one Standard Error \n"),
hjust=0.5) Warning: Ignoring unknown aesthetics: linetype